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SUMMARY 


An order- of-magm'tude analysis of the subsonic three-dimensional, steady 
time-averaged Navier-Stokes equations, for semi-bounded aerodynamic juncture 
geometries, yields the parabolic Navier-Stokes simplification. The numerical 
solution of the resultant pressure Poisson equation is cast into complementary 
and particular parts, yielding an iterative interaction algorithm with an 
exterior three-dimensional potential flow solution. A parabolic transverse 
momentum equation set is constructed, wherein robust enforcement of first- 
order continuity effects is accomplished using a penalty differential con- 
straint concept within a finite element solution algorithm. A Reynolds stress 
constitutive equation, with low turbulence Reynolds number wall functions, is 
employed for closure, using parabolic forms of the two-equation turbulent 
kinetic energy-dissipation equation system. Numerical results document 
accuracy, convergence and utility of the developed finite element algorithm,, 
and the CMC;3DPNS computer code applied to an idealized wing-body juncture 
region. Additional results document accuracy aspects of the algorithm turbu- 
lence closure model. 



INTRODUCTION 


A prime requirement in computational aerodynamics is flow prediction in 
juncture regions formed by the intersection of aerodynamic surfaces, e.g., 
wing-body, wing-winglet, pylon-wing, etc. In most instances of interest, the 
associated flow is three-dimensional, subsonic with variable density, and 
turbulent. The characteristic action of such flows is roll-up of a vortex in 
the plane transverse to the chord coordinate, and mass efflux/influx into the 
boundary layer regions located at some distance from the juncture region. The 
requirement of a numerical algorithm for the juncture flow is to predict the 
associated vortex structure, hence compute a corner drag coefficient, and to 
provide transverse plane velocity boundary conditions for a conventional 
three-dimensional boundary layer analysis of the associated farfield flows. 

The essential key aspects of this problem are illustrated in the geometry 
of the idealized exterior subsonic axial corner, see Figure 1, which has 
received considerable theoretical and experimental attention. Rubin, et al., 
]]l-4] pioneered in formulation and analysis of the three-dimensional laminar 
corner flow problem. Tokuda [5] documents an extension of this analysis, and 
compared his predictions to the experimental data of Zamir and Young £6]. 

Bragg £Zj analyzed the corresponding turbulent flow case, and determined the 
corner distribution of the chordwise Reynolds normal stress component ufuf . 

The salient feature of the turbulent flow case is inducement of a persistent 
axial vorticity component. Various causal mechanisms have been theorized, 
including transverse pressure waves £8] , Reynolds shear stress gradients 
along the corner bisector [g] , and nonisotropy of the Reynolds stress tensor 
£l0] . Quality experimental data for a confined corner flow £ll] , compared to 
documentary results reported herein, indicate the primary mechanism to be noni- 
sotropy of the Reynolds stress tensor "uTuT. 

* vJ 
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Figure 1 Idealized Juncture Region Geometry 


The character in the idealized corner region flow thus appears the result 
of a delicate balance between turbulence phenomena and the induced secondary 
mean flow velocity field. These mechanisms represent a balancing of higher- 
order effects however, as discussed herein, and can be readily dominated by 
flow-field curvature induced vorticity, cf. [12,13], Nevertheless, an adequate 
Reynolds stress closure model is required and has been developed for this 
problem class. The six components of the (symmetric) Reynolds, stress tensor 
are determined using a tensor field constitutive equation formulation which 
requires solution of parabolized forms of the transport equations for turbulent 
kinetic energy (k) and isotropic dissipation function (e). The stress consti- 
tutive equation includes a low turbulence Reynolds number length scale model 
to permit solution of the (k, e) equation system directly through the sublayer 
region adjacent to an aerodynamic surface. Hence, the boundary conditions for 
the k and e solutions are identical vanishing at all aerodynamic surfaces, 

A pressure-velocity formulation is undoubtedly preferred for an algorithm 
to predict turbulent aerodynamic juncture region flows. While definition of a 
transverse plane potential function [14] can automatically satisfy the contin- 
uity equation, the elimination of transverse pressure gradients comes at the 
expense of definition and use of vorticity. The acknowledged weakness of the 
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vorticity formulation is the kinematic boundary condition statement. The 
existence of very large mean flow strain rates at an aerodynamic surface, for 
turbulent flow, serves to further complicate this intrinsic weakness. Con- 
versely, in a physical variable formulation, an algorithm is required developed 
to construct an overall parabolic, i.e., initial -value, elliptic boundary value 
statement for transverse plane phenomena. A careful order of magnitude analysis 
of the transverse plane momentum equations indicates that pressure distributions 
will balance convection and/or turbulence effects to first order, and that 
overall j this balance is of higher order effects than controlled by the con- 
tinuity equation. Since the continuity equation is not parabolic for subsonic 
flow, the construction of a suitable transverse plane equation system is 
required and presented. 

The pressure-velocity formulation is derived and evaluated herein for 
steady turbulent flow prediction in three-dimensional , semi -bounded aerodynamic 
juncture region domains. Persistence of the chordwise component cf the time- 
averaged, mean flow velocity permits an order of magnitude analysis, yielding 
the parabolic approximation to the governing three-dimensional, steady time- 
averaged Navier-Stokes equations. Using the same procedure for components of 
the Reynolds stress tensor, the balancing of lowest order terms in the two- 
transverse momentum equations yields a pressure Poisson equation. An algo- 
rithm for this equation is derived in terms of complementary and particular 
solution fields. The complementary solution is determined using boundary 
conditions obtained from an exterior potential flow solution. The particular 
solution refines this pressure determination by accounting for the Reynolds 
stress and transverse velocity distributions. The particular solution is 
enforced in a retarded manner in the chordwise momentum equation, to update 
the three-dimensional pressure field, yielding an iterative-interaction algo- 
rithm with, the three-dimensional exterior potential flow solution. Algorithm 
convergence occurs when this composite pressure solution becomes stationary. 

As a consequence of the ordering analysis, the number of dependent variables 
requiring solution exceeds the available equations. Therefore, using finite 
element penalty function concepts in constrained extremization, cf. [15] , 
a transverse momentum equation solution statement is constructed wherein the 
first-order effects of the nonparabolic continuity equation are enforced as 
a differential constraint. 
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SYMBOLS 


a 

A 

b 

B 

C 

E 

f 

F 

9 

G 

h 

H 

i 

j 

J 

k 

i 

L 

M 

n 

N 

P 

s,S 


"rj 

u 

y* 

a 

9R 

6 

6Q 


boundary condition coefficient 

constant 

constant 

finite element two-dimensional hypermatrix prefix 
turbulence model coefficient; initial-value matrix 
energy norm 

function of known argument 

finite element matrix; discretized equation system 

function of known argument 

finite element matrix prefix 

metric coefficient 

stagnation enthalpy; Hilbert space 

index 

index 

Jacobian matrix 

turbulence kinetic energy; finite element basis degree 

summation index; differential operator 

differential operator 

number of finite elements spanning R^ 

unit normal vector; dimension of space 

finite element cardinal basis; discrete index 

pressure; iteration index 

heat flux vector; generalized dependent variable 
generalized semi-discrete dependent variable 
spatial domain of differential operator 
source term 

finite element assembly operator 
velocity vector 

Reynolds kinematic stress tensor 
convection matrix 
Cartesian coordinate system 
shear velocity Reynolds number 
Lagrange multiplier set 
partial derivative operator 
boundary of solution domain R^ 

Kronecker delta; parameter 
iteration vector 
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A mesh measure; increment 

e isotropic dissipation function; parameter 

n- curvilinear coordinate system 

K heat conductivity coefficient 

V kinematic viscosity 

p density 

0. . Stokes stress tensor 

* J 

Z summation 

(j) constraint dependent variable 

0) sublayer damping function 

£2 solution domain 


Superscripts 


e 

h 

0 

P 

T 

t 


finite element reference 
solution approximation 
initial condition reference 
iteration index 
matrix transpose 
turbulent reference 
ordinary derivative 


Subscri pts 


e 

i,j,k, 5 , 

j 

0 


finite element reference 
tensor indices 
time step index 
reference state 


Notation 

{ } 

[ ] 
u 
n 


column matrix 
square matrix 
union 

intersection 
belongs to 
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PROBLEM DESCRIPTION 


Parabolic Navier-Stokes Equations 

The three-dimensional parabolic Navier-Stokes (3DPNS) equations are a 
simplification of the steady, three-dimensional time-averaged Navier-Stokes 
equations. In Cartesian tensor notation, and employing superscript tilde and. 
bar to denote mass-weighted and conventional time-averaging respectively, 
the conservative equation form for a variable density, heat conducting fluid 
is 


- “ij] = ° 

\J 

U(pR) - |j7[pS0j - u,;,j * 5fP=J - * qj » 0 


* qqfqj 


la + pp = 0 


L(p£) = 


- — pu.e + C 




(1) 

( 2 ) 

(3) 


(4) 


+C^^ = 0 
e k 


(5) 


In equations 1-3, p is density, u. is the mean velocity vector, p is pressure, 
is the Kronecker delta, and H is stagnation enthalpy. The Stokes stress 
tensor a., and heat flux vector q. are defined as, 

* J J 


°ij “ - 




q 


j 



( 6 ) 

(7) 
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where Re is the reference Reynolds number, Re = U L/v , and -puiu' is the 

00 OO 1 J 

Reynolds stress tensor. In equations 6-7, v and k are fluid kinematic viscos- 
ity and heat conductivity respectively, and E.. is the mean flow strain rate 

* J 

tensor 



( 8 ) 


Equations 4-5 are the transport equations for turbulent kinetic energy and 
isotropic dissipation function, as obtained using the closure model of Launder, 
Reece and Rodi [17} for the pressure-strain and triple correlations, and 

k = (9) 


= ^^1 r 

^ ■ 3 3Xj jk 


( 10 ) 


The various coefficients C“ are model constants, cf. [18]. 

3 1- -I 

The parabolic Navier-Stokes equation set is derived from equations 1-5 by 
assuming the ratio of transverse mean velocity components to chordwise 
component is less than unity, and by further assuming that: 

1 . the chordwise velocity component suffers no reversal, 

2 . diffusive transport processes in the chordwise direction 
are higher-order, hence negligible, and 

3. the overall elliptic character of the parent three- 
dimensional Navier-Stokes equation is enforcable through 
construction of a suitable pressure field with exterior 
flow boundary conditions. 

Assume the x^ (curvilinear) coordinate direction parallel to the chordwise 
mean flow direction, with scalar velocity component Uj of order unity, i.e. 
0(1). Further assume 0 (u 2 ) ~ 0 ( 6 ) ~ 0 ( 03 ), and that 0 ( 6 ) < 0(1). As occurs 
with boundary layer theory, the continuity equation confirms that chordwise 
variation in Ui is of the order equal to appropriate transverse variation of U 2 
and U 3 ; hence, for = 0 ( 1 ), =» 0 ( 6 "M = • 
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Determination of the relative order of terms in the momentum equation 2 

is straightforward. For the u^ equation, since O(puTuT) must be 0(6), the 

term -r— (pufuf) is higher order and can be discarded. The assumption that 
^ ^ 1 9 ^ 1 
chordwise diffusion is negligible infers that (En) = 0, hence 0(Re~ ) £ 

0(6). Therefoire, the terms in 5i2 and 5 i 3 involving U 2 and D 3 i.e., 

■^^- 1 , are both 0 ( 6 ) or smaller and can be neglected. 

3Xij 

Deletion of these terms is fundamental to the parabolic approximation, 
since their elimination removes the elliptic boundary value character in the 
chordwise flow direction. The existence of -r— (pUiUi) instills an initial 

oX 2 

value character in the resultant equation, hence, permits marching the solution 
for Ui in the chordwise direction. The desired 3DPNS form, denoted lP(*)> is 
therefore. 


3 Xo 


and ’ 


9Xi 


3Xa 


lP(p-uO . (SufuO . 4^ + 


3X2 


pUjUf - O12 


+ 




(11) 


which is thoroughly familiar. As a final note, should x- correspond to a 

J 

curvilinear coordinate description, the derivatives expressed in equation 11 
are interpreted as covariant derivatives. The 3DPNS form of the energy 
equation (3), similarly constructed, is 




pHUc 


“ U .a . 

1 U 


+ pH'u: 


- U'o' 

1 u 


= 0 


( 12 ) 


Equation 12 introduces the 3DPNS limited index summation convention 
1 £ (i»j) £ 3 and 2 £ £ ^ 3. 

In agreement with boundary layer concepts, the order of pressure variation 
in the transverse plane is assumed controlled by the lowest order terms appear- 
ing in equation 2 written on U 2 and Ug. Each transverse derivative of pU£Uj 
and puIuT is 0(1), while all other terms are 0(6) and higher. Thus, for a 
conventional two-dimensional boundary layer flow, for example, 


lP (pUz) 


3X2 


P 


PU 2 U 2 


= 0 


(13) 
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The solution is trivial; p differs from the inviscid flow edge pressure by a 
constant, equal to a fraction of the free-stream turbulence (k) level, and is 
distributed through the boundary layer in proportion to FuJuJ* The initial 
value character for pressure, as exhibited by the 3DPNS first order approxi- 
mation to equation 2,for D2 and U3 , is recast into a more tractable form by 
taking the divergence. Retaining the higher-order convection and diffusion 
terms for generality, the consistent 3DPNS form for both transverse momentum 
equations is 


•-(P) = Irf ■^1 


3X2 3X.3X„ 

Z J £ 


pu„u,. + p uru; - a„ 

_ J!- J £ J £jj 


= 0 


(14) 


Equation 14 defines an elliptic boundary value problem for determination 
of pressure distribution in the transverse plane. The pressure field that 
satisfies this quasi-linear Poisson equation consists of complementary and 
particular solutions, i.e.. 


P(x,) . p^(x,) + Pp(x,) 


(15) 


The complementary solution satisfies the homogeneous form of equation 14, i.e.. 


lP(p^) 


3X^ 


= 0 


(16) 


The Dirichlet boundary condition for equation 16 is P^,(xi, x^) = p(xi, x^), on 
the intersection of the 3DPNS domain with the exterior potential flow domain . 
Elsewhere, the boundary condition for p^ is homogeneous Neumann. 

The particular pressure is any solution to equation 14 subject to homogeneous 
Dirichlet boundary conditions on boundary segments where p^ is known. Elsewhere, 
the nonhomogeneous Neumann constraint is provided by the inner product of the 
3DPNS form of equation 2, written on u^^, with the local outward pointing unit 
normal ri^^. 


£(p ) E LP(pu.)-n^ = 


9x, 


Pp ^ P V£ 


3Xi 


P^k"£ 


= 0 


(17) 


10 



Repeated indices in equation 17 are not summed, and 2 k ^ 3 for k 
Hence, equation 17 is the generalization of the boundary layer form, equation 13. 
Following determination of the order of terms in the Reynolds stress tensor in 
the next section, the nonhomogeneous terms in equation 17 vSnish to lowest 
order on an aerodynamic surface. Hence, thereupon p_(xi, x ) is a constant 
which is zero. Elsewhere on the 3DPNS domain boundary, equation 17 yields the 
appropriate boundary condition for equation 14. 

Reynolds Stress Tensor Closure 

A closure expression for the kinematic Reynolds stress tensor -lI^uT , 
appearing in equations 4-5, is required to complete the 3DPNS order of 
magnitude analysis. The necessary insight is provided by construction of a 
tensor field strain-rate constitutive equation, the existence of which is 
assured at "sufficient" distance from boundaries in space and time [19]. 

Using lower-dimensional order of magnitude analyses and invariance, the 
three lead terms of the five term expansion of the kinematic form, appropriate 
for 3DPNS analyses , are [ 20 ], 


-u.u. 
1 J 


= -ka 


TJ 


k^ 


+ C 




(18) 


E.. is the symmetric mean flow strain-rate tensor, equation 8, and k and e 

1 vi 

are turbulence parameters defined in equations 9-10. 

Equation 18 results from re-expression of the triple correlations, within 
the Reynolds stress transport equation, using the model of Launder, Reece and 
Rodi [17] , and is the tensor generalization of the original analysis by 
Gessner and Emery [lO] . In equation 18, a. , is a diagonal tensor in the prin- 

• J 

cipal coordinates, defined as 


‘ij ~ 3k ^'^k‘^k^®i'^ij 


(19) 


The a^ are coefficients admitting anisotropy, where a^ = Ci, and a 2 = C 3 = as. 
The C” are defined [17] as 
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22(Coi - 1) - 6(4 Co2 - 5) 

- 33(Coi - 2 Co2) 

4(3Co2 - 1) 

= lUCoi - 2 Co2) 

22(Coi - 1) - 12(3Co2 - 1) 

- 33(Coi - 2 Co2) 

44Co2 - 22CoiCo 2 - 128Co2 - 36Co2 + TO 
= 165(Coi - 2 Co2)^ 


( 20 ) 


In equation 20, Coi and Cq 2 are "universal" empirical constants; suggested 
values are Cqi s 2.8 and Cq 2 = 0.45, [18]. 


The order of terms in equation 18 can be estimated for the standard 

values C = {0.94, 0.067, 0.56, 0.068}. For significance in equation 2, recall 
a 

OCu'u') = 0(6). For a two-dimensional flow, and for i = 1 and j = 2 , equation 


18 yields the familiar form 


-UfU2 = Cl 


k2 


8U] 

9X- 


( 21 ) 


Hence, 0 (C 4 k 2 /e) = 0(62). Further for i = 1 = j, and neglecting the second 
two terms, 0 (k) = 0 ( 6 ), hence 0 ( 04 / 0 ) ~ 0 (i)* To proceed further requires 
an estimate of the magnitude of 0(k). For a steady, subsonic, turbulent aero- 
dynamic flow, away from the influence of the wall or freestream, the magnitude 
of the fluctuating component of velocity will probably not exceed about 10 % of 
the steady component. Hence, equation 9 in nondimensional form yields the 
estimate 0(k) £0(10" ). Taking the maximum yields 0(6) £0(10” ); evaluating 
the fourth term of equation 20 and comparing yields 0 ( 64 ) = 0 ( 6 '^), hence 
0(e) = 0(6'^). Thus, 0 ( 0264 -^ ) 2 0(62). Therefore, in rectangular Cartesian 
coordinates, and retaining terms of the first two orders of significance, the 
six components of the kinematic Reynolds stress tensor for a 3DPNS analysis are 
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r- 


0 ( 6 ) 






0 ( 6 ") 




UiUf = Cik - C2C4— T 

£ 


3Ui " ^ 

’3U1 

2 

3 X 2 

[ 3 X 3 J 



Ua'Uj = Csk - 


b 3 n^ri_~l 2 

CaC^pl 


= Cak - 


UfU2 = 


- C, 


f auil ' 

[3X2j 

Psui j ' 

[axaj 

k" r^uil 


k" 

C2Ci»-^ 


U1U3 = 


- Cu-^ 


U2U3 = 


C2C4- 


- 2C, 


- 2C, 


- ZCu 


- 


TaOi' 
e _9Xi_ 



'^r3U3~! 

-9X3J 

9Ui r3U2 ^ aual 

3x31,3x3 9X2J 


+ 2 - 


3x2 


3Xi 3X2 


3Ui 


3Ui 

f3U2 + 3U3l 

_3X3j 

U 0 U L. y 

e 

3 X 2 

, 3 X 3 3 X 2 ; 


+ 2-^f-^ + 

3 X 3 ^3Xi 3XaJ 


3Ui 3Ui 


3U2 _|_ 9^3 

_3X2 3X3_ 

- 

3 X 3 3X2_ 


( 22 ) 


Two conclusions regarding equation 22 should be noted. The terms which 
would provide an elliptic boundary value definition in the X2« x^ plane, for 
direct integration of equations 2, for U2 and Q3, cf. [21], are indeed 0(6^), 
in agreement with the ordering arguments leading to equation 14. Secondly, the 
0(6) term In -ufiJJ In equation 22 vanishes on an aerodynamic surface, hence 
thereupon pp is a constant. With this development, the order of terms in equa- 
tions 4-5 can be determined, yielding the appropriate 3DPNS approximation as. 


'■'(k) - - 5)|^' 


* 5 ^ - 0 


(23) 
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(24) 


lP(£) = 




1 L 


C iiTJrUT 1^1 

eel a 3Xj^ 


^e kl^ 


+ C^ p I- =0 


recalling the 3DPNS limited summation convention, 1 i £ 3, 2 ^ z ^3. 

Definition of boundary conditions for equations 23-24 requires addressing 
the issue of what constitutes "sufficient" distance for validity of equation 
18. For two-dimensional flows, one approach is to employ similarity arguments 
to assign values to k and e at some distance from the wall, e.g. 10 < y'*’ < 50, 
where y*" £ ''s a turbulence Reynolds number based on wall shear velocity 

u^ = A^/p , see [l6j. Extension of this concept to three-dimensional flows 
is questionable, but has been attempted l_2l] . A second alternative [22] 
suggests modifying the "constants" C“ appearing in equations 23-24, and 

P 

integrating directly through the low turbulence wall region with k £ 0 £ e as 
boundary conditions. 

The alternative approach of £23] is employed for the juncture region 
analysis. The "constants" of the Reynolds stress constitutive equation 
18 are modified to account for low turbulence levels in the sublayer region. 
Equation 21 defines the conventional turbulent "eddy viscosity"' h C 4 k^/e. 
Using dimensional analysis, is the product of a scale velocity and a scale 
length; typically, for a turbulence kinetic energy model, 

= k"£^ (25) 


Comparison with equation 21 yields the familiar relationship 


£ 


d 



(26) 


Recalling the van Driest damping function u, defined to control evolution of 
the Prandtl mixing length scale £l6] , equation 26 multiplied by oj yields 
equation 21 in the form 


-u{u£ 



3ui 

3X2 


(27) 
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The conventional form [16} for u is modified,, for variable length scale damping, 
as 

OJ = [j - exp(-byVA'^)j (28) 

where A"*^ » 26, and b « 2.0 based upon results of numerical studies £23, 24]. 
Therefore, premultiplying each of the coefficients in equation 20 by to 
produces the required sublayer modification for equation 18. Furthermore, 

C| in equation 24 is also multiplied by w. The aerodynamic surface boundary 
conditions for equations 23-24 are then k = 0 e e. 

Differential Eq_uation System Closure 

Development of the lowest order parabolic Navier-Stokes differential 
equation system, as a subset of the steady, time-averaged Navier-Stokes system 
is complete. This 3DPNS system, equations 11, 12, 14, 18, 23 and 24, numbers 
one less than the number of dependent variables defined in equations 1 - 10 . 
Therefore, at least one equation governing 0(6) phenomena must be included to 
close the system. Since the 3DPNS momentum equation 11 is written on u^ only, 
both components of u^ = { 02 . 03 } are required determined subject to the con- 
straint of continuity, equation 1. The finite element algorithm accomplishes 
this by "penalizing" the solution of the 0(6) 3DPNS approximation to the momen- 
tum equations 2 , written on both components of u^, by the continuity equation 
(error). Retaining the first two orders of terms, the 3DPNS form for the 
transverse momentum equations is. 




(29) 


9X 








which introduces the additional 3DPNS limited index 2 £ k _< 3. The middle two 
terms in the second bracket are 0 ( 6 ), while the remaining terms are all 0 ( 6 ^) 
or smaller. Equation 29 exhibits elliptic boundary value character in the 
Xg. Xg plane, retaining the terms of 0(6^) in the Reynolds stress tensor, 
equation 18, and contains the initial-value term permitting chordwise marching. 
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Since equation 29 represents two additional scalar equations, an auxiliary 
dependent variable is required defined. The theoretical concept, borrowed from 
the variational calculus, is to define a suitable measure of the continuity 
equation (error), which is then applied as a differential constraint on solu- 
tion of the transverse momentum equation 29. This constraint measure must span 
the transverse plane R^, and must vanish as the continuity equation 1 becomes 
satisfied. Based on computational experience [23], an appropriate dependent 
variable is the harmonic function defined as the solution to the Poisson 

equation 



subject to homogeneous Neumann boundary conditions on portions of the domain 
boundary 8R, and setting <j> = 0 at one location at least on aR. Equation 30 
becomes homogeneous, as the continuity equation becomes satisfied, and the 
solution 4>(x^) becomes null as a consequence of the boundary condition 
specifications . 

Grid Stretching Transformation 

An elementary grid stretching coordinate transformation is of potential 
use for the general problem class, and is consistent with the ordering simpli- 
fications yielding the 3DPNS equation system. The transformation n^. = n^-(Xj ) 
that normalizes transverse spans with boundaries f^^- , 2 £ 5 , _< 3, 1 £ i £ 2, is 


{n,} 


X2 ~ f2l(Xl) 

[f22(Xl) - f2l(Xl)]/f2 

Xa ~ faiCXj) 

[f32(Xl) - f3l(Xl)]/f3 


(31) 


The f^^-(xi) are piecewise continuous segments defining the transverse plane 
boundary 3R of R^, and the f are normalizing coefficients. Using the chain 

X/ 

rule, differentiation on introduces additional derivatives on In 

particular. 
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- [hz2 + nahas] - [fi 32 + nahas] 


d ^ 9 

3 Xi dVii 


3 = \y 

( 32 ) 


The functions 1 ^ 5. defined as 


<v,) =- 


h„f:" (f^2 - 


V -£l-£ 


(33) 


The superscript prime denotes the (ordinary) derivative with respect to 
and the coordinate system is fixed in the transform space. 
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FINITE ELEMENT SOLUTION ALGORITHM 


The consistently ordered 3DPNS equation system has been constructed for the 
dependent variable set q-(x.) e {q} = {p, Uj, U 2 , U 3 , H, p, k,e, u^ui, 

An equation of state p = p(p, H) closes the system. Equations 11, 12, 23, 24, 
and 29 of the 3DPNS equation set contain the initial -value term that facili- 
tates solution marching in the chordwise direction. Equations 14 and 30 are 
elliptic boundary value descriptions with parametric initial -value dependence, 
while equation 18 is a local constitutive definition. Equation 1 becomes 
recast as the differential constraint using equation 30. 

The general form of 3DPNS system description is 



For equation 34, g • and s. are specified non-linear functions of their argu- 
ments, as determined by the index j. The three-dimensional partial differen- 
tial equation 34 is defined on the Euclidian space spanned by the x^. (,n^. ) 
coordinate system. The solution domain o. is defined as the product of R^ 
and Xj, for all elements of x^ belonging to the open interval measured from 
Xi( 0 ), i.e., 

Q E r2 X Xi = {(Xj^,xi): x^eR^ and Xie[xi (0) ,Xi) } 

the boundary of the solution domain is the product of the boundary 9R of R^ 
and Xi, i.e., dn = 3R x x^. Thereupon, a differential constraint is applied 
of the form 

^(Qj) = aiq^ + a2 q^.n^. + as = 0 (35) 

In equation 35, the a^. are specified coefficients and h^. is the outwards 
pointing unit normal vector. Finally, an initial distribution for the 
appropriate members of q. on f 2 o = R^ x Xi( 0 ) is required. 

J 



For the finite element numerical solution algorithm of equations 34-35, the 

approximation q!-(x , x^) to the (unknown) exact solution q-(x , x^) is con- 

structed from members of a finite-dimensional subspace of H^(n), the Hilbert 

0 

space of all functions possessing square integrable first derivatives and 
satisfying the boundary condition 35. While extremely flexible in theory, the 
practice for the 3DPNS equation system is to employ linear polynomials, 
defined on disjoint interior triangular-shaped subdomains R|, the union of 
which forms the discretization of R^. Hence, the finite element approximation 
is 

qj(x^,X:) = q 5 (Xj,xO = q-(Xj.Xi) ( 37 ) 

using the elemental construction 

q 5 (x^,xi) i (Ni(Xj)}'’’{QJ(xi))g ( 33 , 

L 

In equations 37-38, j is a free index denoting members of {q^}, and subscript 
or superscript e denotes pertaining to the e— finite element, = R^ x xi. 

T G G 

The elements of the row matrix {Ni(x„)} are linear polynomials on x , 2 < 

a z — 

z < 3, £28], and elements of {QJ}„ are the values of q. at the nodes of R^ . 

— e j e 

The functional requirement of any numerical solution algorithm is to render 
the error in q^ minimum in some norm. The finite element algorithm requires 

J n h h 

the error in equations 34 and 35, i.e., L^(q^) and to be orthogonal 

h J J 

to the space {Ni(x )} employed to define q-. In addition, the discrete approx- 

D -h ^ J 

imation L'^(p ) to the continuity equation 1 must be enforced as a differential 
constraint. Identifying the (Lagrange) multiplier set 3 .j , these linearly 
independent constraints are combined to yield the finite element solution algo- 
rithm theoretical statement. 

/R 2 {N^}LP(q 5 )dx+ 8 i/g|^{N^H(q 5 )dx + l 2 ./p 2 V{N^}LP(p^d^ = 

Equation 39 represents a system of ordinary differential equations, written 
on the chordwise coordinate x^, of the ■form 

£C]{QJ}" + £U]{QJ} + [ 6 LJKQL} + {SJ} = {0} (40) 
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A one-to-one correspondence of terms in equations 40 and 34 is inferred, as 
augmented for the various additional terms introduced through 3^ 0 in 

equation 39. The integration algorithm for equation 40 is the trapezoidal 
rule; hence. 


{FJ} = - {QJ}j = {0} 


(41) 


defines a system of nonlinear algebraic equations for determination of the 
elements of {QJ(x2)}. A Newton iteration algorithm is employed for solution 
of equation 41 as 




(42) 


The dependent variable in equation 42 is the iteration vector, related to 
the solution in the conventional manner, 

P+1 p p+1 

{QJ} = {QJ} + {6QJ} (43) 

j+1 j+1 j+1 


v{N|^}LP(p'^)dx 


The algorithmic embodiment of the differential constraint concept 
employs a sequential summation into the column matrix denoted I2 
in equation 39. A numerically determined optimum expansion coefficient is 
32 = Axi{j,k}, where j and k are unit vectors parallel to x^. In the trans- 
verse plane momentum equations, this term corresponds to a load (column) 
matrix, say {G2PHI} for U2. Letting denote the nodal solution for 


(j)' (x ) at iteration step i then 

Xf 


{G2PHI}j^^ E j.V^PHI}]^^ (44) 

where denotes the integral of the discrete gradient operator on the mesh 
of measure h. This contribution is then added to the sum of the previous p 
evaluations for cj) , to construct the action of the differential constraint 
term for step (AXi)j_|_-|, iteration p+1, i.e., 

p+1 p i 

{GU2}.,, E Z {G2PHI} (45) 

^ ' i=l j+1 
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Hence, each successive determination of {PHI} corrects the action of all 

j+1 

previous solution iterates, such that {PHI}P_|_-|->-{e} as p increases without 

bound, where \z\ > 0 is an acceptable discrete level of computed zero. This 

procedure thus admits, in the limit, the exact continuity preserving solution 
for equation 30, i.e., ^ ~ 0 everywhere. Additional discussion on details of 
algorithmic constructions are given in the Appendix. 


THEORETICAL ANALYSIS, ACCURACY AND CONVERGENCE 


The 3DPNS equation system contains as a subset the two-dimensional 
boundary layer equations for laminar or turbulent flow. For these elementary 
systems, a finite difference truncation error analysis confirms the linear 
basis finite element formulation is spatially second-order accurate. Of 
course, the trapezoidal rule employed for chordwise marching is also second- 
order accurate. A formal analysis of convergence in Sobolev norms, [26] for a 
scalar linear parabolic equation,- predicts the error in the semi-discrete 

L. 

linear basis finite element approximation q” satisfies the inequality. 


[e^(nAXi ) ,e^(nAX]^ ) 


< CiA! 


|q(nAxi) I 1^ + C2 Ax2 


(46) 


where C^ is a constant independent of A^, the measure of the largest finite 
element on R . Furthermore, C 2 is a constant independent of a^, axi is the 
space-marching step, and ||Qq|P is the "energy" in the initial data. Hence, 
equation 46 confirms the solution error is bounded by a constant times a term 
of order a|, i.e., second-order accurate. Furthermore, from the fundamental 
theorem j[27] , the semi-discrete approximate solution converges in energy, i.e., 

ECe*^, ^ 0 as Ag 0 (47) 

The strongly nonlinear 3DPNS (and turbulent boundary layer) differential 
equation systems are significant departures from the elementary equations 
considered in ^.26, 27] . For example, the energy norm E(-,-) is evaluated as 
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( 48 ) 



The "effective" turbulent diffusion coefficient is a strongly 

non-linear function of with both k and e exhibiting nearly singular 
behavior in the sub-layer region immediately adjacent to an aerodynamic 
surface. Nevertheless, the results of closely controlled numerical experi- 
ments [28, 29 J have predicted extension of the linear theory, for the boundary 
layer equations, as well as providing exact comparisons between the k = 1 


finite element algorithm, and the equal complexity (and familiar) Crank- 
Ni col son finite difference algorithm. 


The fundamental theoretical aspect of critical importance is quantization 
of performance of the penalty continuity constraint on the transverse momentum 
equation solutions. The classical concept of the penalty construction for a 
linear elliptic statement [26^ defines the parameter ^2 the (norm of 

the) penalty term approaches zero. For the nonlinear parabolic Navier- 
Stokes system, and following extensive numerical experimentation, 

“ AXi[j+kj was determined preferable in optimizing the number of iterations/ 
step to convergence. Using the outlined accumulation procedure, equations 
44-45, the penalty algorithm yields satisfaction of the continuity equation 
in energy to the order of E((j)^,(f)'^ ~ 0(10"®) for laminar boundary layer flow, 
E(<f)*^, (j)*^) = 0 ( 10 "^) for turbulent boundary layer flow, and E((|)*^, 4>*^) = 0(10 ^) 
for 3DPNS solutions. 
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PROBLEM ANALYSES 

Idealized Wing-Body Juncture Flow 

The documentary test case configuration, for the developed 3DPNS algorithm 
and the CMC:3DPNS computer program, is turbulent flow in the juncture region 
formed by the right intersection of two 10-percent thick parabolic arcs with 
coincident leading edge. A complete discussion of the CMC:3DPNS program 
input and data deck preparation procedures is presented in Volume II of this 
report. The computer program documentation manual for the CMC;3DPNS code is 
published as Volume III. 

Figure 2a illustrates the essential geometry of the parabolic arc juncture 
specification. The complementary pressure (p^) boundary conditions for the 
3DPNS solution were obtained using the Hess potential flow computer code [30], 
for M = 0.08 (U = 30 m/s = 100 f/s) and zero angle of attack. Figure 2b) 

OO CO 

summarizes the resultant spanwise distributions of P_(x ) at chord stations 
Xj/C = 0.01, 0.085, and 0.46. By symmetry, these pressure boundary conditions 
are also appropriate at x^/C = 0.54, 0.915, and 0.99. Therefore, a progres- 
sively decreasing favorable x^ pressure gradient exists to mid-chord, and 
thereafter turns progressively adverse. The strongest gradients are confined 
to the immediate vicinity of the corner. 

The 3DPNS solution domain was defined to span 0 ^ x^/C ^ 0.1 and of height 
X|^/C » 0.01. Figure 3 illustrates a nonuniform discretization of the trans- 
verse plane as the union of triangles (with most diagonals omitted for 
clarity). The domain boundary 3R is the union of straight line segments A-F, 
upon which boundary condition specifications are required for the dependent 
variable set q^(x ,xi), see Table 1. The Reynolds number is Re/C = 0.6 x lO^/m, 
and the flow is assumed isoenergetic, hence H(x^) = constant. The initial 
conditions for u°(x ), at the nodes of R^, are established using Cole's law to 

Xf 

interpolate a turbulent boundary layer profile onto node "columns," Figure 3, 
with matching of the free-stream level of p (x ). The transverse plane velo- 
city Uo(x ) is defined identically zero, until eight 3DPNS steps are completed, 
to permit computation of a reasonable chordwise derivative of pui, e.g., 
{RHOUl}^. The initial distributions for k°(x ) and e°(x ) are computed using 

36 36 

boundary layer mixing length concepts, as discussed in detail in [23]. Each 
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Span Coordinate X2/C 






TABLE 1 


Boundary Condition Specifications for 
3DPNS Analysis of Aerodynamic Juncture Region Flow 


Dependent 

Variable 

Equation 

No. 

Boundary 

Seqment 

Boundary Condition 
Statement, Equation 35 

Comments 





ai 

a2 

^3 




11 

ABC 

1 

0 

0 

No-slip Surface 




C.DEFA 

0 

1 

0 

Vanishing Derivative 



29 

ABC 

1 

0 

n 

V.' 





CDEFA 

0 

1 

0 



k 

23 

ABC 

1 

0 

0 





CDEFA 

0 

1 

0 



e 

24 

ABC 

1 

0 

0 





CDEFA 

0 

1 

0 



p 

16 

DEF 

1 

0 

f) 

Potential Constraint 



J- KJ 

FABCD 

0 

1 

0 



Pp 

14 

ABC 

1 

0 

0 



! 


CDEFA 

0 

1 

0 



4> 

30 

DEF 

FABCD 

1 

0 

0 

1 

0 

0 

1 




of these initialization procedures is a CMC;3DPNS input specification, 
discussed in detail in Volume II. 

The standard test case specification is 3DPNS solution of the parabolic 
arc juncture region turbulent flow on the interval 0.01 £ x^/C _< 0.60, using 
the 10 X 19 node nonuniform discretization shown in Figure 3. The corres- 
ponding number of triangular finite element domains is M = (1 0-1 )X(1 9-1 )X2 + 1 
= 325. On this grid, the 3DPNS algorithm takes 170 chordwise steps and 
averages 4.3 Newton iterations/step for convergence set at e = 0.0003. On 
the NASA Langley CDC Cyber/203 computer, with no vectorization of the scalar 
CMC:3DPNS code, this execution requires 475 seconds of CPU time. For com- 
parison, the same execution on an IBM 370/3031 computer requires approximately 
6000 seconds of CPU time. The central memory requirement for both executions 
is 150,000 words. During the 3DPNS solution, the distribution of the particu- 
lar pressure solution Pp(x^,Xi) is written on an output file at each chord- 
wise station x^ for which the p^(x , x^) boundary condition is specified. 

For the second and sequential 3DPNS solutions, the solution of equation 16 
for Pc(x^,Xj) is algebraically summed with the stored distributions Pp(x^, 
x^), see equation 15, and this sum employed for the chordwise pressure 
gradient distribution for the u^ momentum equation solution, see equation 11. 
The u^ momentum equations are solved using the current computed Pp distribu- 
tions. The composite pressure field p(x.j)/Po, equation 15, converges to five 
significant digits following three 3DPNS algorithm solutions, for the standard 
test juncture region geometry. The nominal level of Pp/Pg is 10'^; a repre- 
sentative extremum difference between the second and third 3DPNS solutions is 
APp/Pp ~ 10“^ at Xi/C = 0.17. 

The 3DPNS solution computes and outputs the distribution of q.(x ,xi) at 

vJ 

select chordwise stations Xj. Figures 4-5 illustrate qualitatively the third 
interaction 3DPNS solution for transverse velocity for the parabolic arc 
juncture. Figure 4 graphs the complete transverse plane velocity distribu- 
tion Uj;,(xj^) at Xi/C = 0.31 on the M = 325 grid. The solution is exactly 
mirror symmetric, a generated vortex pair in the corner is just distinguish- 
able, and the solution computes a spanwise efflux of u^ from the 3DPNS domain. 
Figure 5 summarizes the M = 325 grid solution evolution of the lower surface 
distribution of Uj^(x^) on 0.021 _< x /C ^ 0.7. (The truncated upper portions 
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Figure 4. 3DPNS Solution for Transverse Plane Velocity u^ Distribution, 
Parabolic Arc Juncture Region, xi/C= 0.40, M = 325 , Turbulent Flow. 
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0.000 0.0,17 0.034 0,. 051 0.06S C.CS5 0.102 


a) Xi/C = 0.021, u'J/u^ = 0.167 



0.000 0.017 0.034 0.051 0.068 0.085 0.102 


b) Xi/c = 0.047, u'J/Ui = 0.072 



c) Xi/C = 0.081, = 0.058 



0.000 0.017 0.034 0.05.1 0.06.8 0.085 0.102 


d) Xi/C = 0.173, u'J/iji = 0.112 

Figure 5. 3DPNS Solution for Transverse Plane Velocity u^^ 
Distributions, Parabolic Arc Juncture Region, Turbulent Flow. 
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0.000 0.01?' 0.034 0.051 0.068 0.085 0.102 


e) .Xi/C = 0.349, u"*/ai = 0.114 

^ / \ \ \ \ N '--'X 


^ S N N N, ri. 

•n N. > — >• >■ 



0.000 0.017 0.034 0.051 0.068 0.085 0.i02 


g) Xi/C = 0.631, u'J/Oi = 0.098 



0.000 0.017 0.034 0.051 0.068 0.085 , C . I 02 

h) xi/C = 0.704, u'J/ui = 0.097 


Figure 5. Concluded. 
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are exactly mirror symmetric.) A massive influx into the corner region is pre- 
dicted at Xi/C = 0.021, with the extremum scalar component u'^/u = 0.167, i.e., 
equal to 17% of the local free-stream value of Ui. This is the direct result of 
the associated large favorable pressure gradient for u^, coupled with the fact 
that the mass conservation algorithm has just become initialized. By Xi/C = 
0.047, the extremum scalar component is = 0.072, and the juncture blockage 

is inducing a large spanwise efflux from in the lower reaches of the boundary 
layer. A corner axial vortex (pair) is just visible at X 2 /C = 0.081, where the 
minimum level of u is predicted. This general velocity distribution persists 

A/ 

to well beyond mid-chord, with “ 10% throughout. The corner vortex pair 

becomes fully developed, with the spanwise efflux filling the entire boundary 
layer. Past mid-chord, the free-stream velocity derivative changes sign, with 
a concurrent cessation of influx from the potential region into the corner 
indicated by x^/C = 0.63. 

A repeating of this solution for laminar flow provides an additional quali- 
tative accuracy assessment. Figure 6 compares the 3DPNS transverse plane 
velocity distributions at xi/C = 0.46, for laminar and turbulent flow solutions. 
In comparison, the corner vortex pair is slightly larger for laminar flow, the 
extremum component u'J'/ui = 0.06 is 40% smaller in magnitude, and a reversed 
spanwise flow is predicted in the lowest reach of the farfield boundary layer. 
Figure 6a). Shafir and Rubin [31] predict theoretically this lower reach 
reversal, for laminar-turbulent boundary layer transition, see Figure 7. 
Furthermore, the flow-fields in Figure 6 are in qualitative agreement with the 
composite corner layer/asymptotic boundary layer solution reported by Rubin 
and Grossman [3] . 

Figure 8 is a composite summary of pertinent transverse plane isoclines 
of Ui and components of u^u^ at xi/C = 0.46, as predicted by the third inter- 
action 3DPNS solution. The Ui solution exhibits the intrinsic symmetry with 
a modest relative displacement directly adjacent to the corner, the result of 
the axial vortex pumping low momentum fluid into the corner and out parallel 
to the diagonal. Figure 6(b). The plot of TITuT is also symmetric, and a local 
extrema exists in the corner due to the axial vortex pumping of the wall layer 
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a) Laminar Flow, u^/ui = 0.06. 



Figure 6. 3DPNS Solution Transverse Plane Velocity u^^ 
Distribution, Juncture Region Flow, Xi/C = 0.46. 



Figure 7. Boundary Layer Laminar/Turbulent Transition Solutions 

In Juncture Region Farfield, from Shafir and Rubin [31] 
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flow into the corner. These comments are valid for the distribution of ulujs 
except that this normal stress exhibits a (very) modest nonsynmetry due to the 
0(6) terms involving Ux in equation 22. The ufuj shear stress distribution is 
highly honsymmetric with the extremum level penetrating only half the span 
distance to the corner. On the upper-half domain, these levels are very small 
since Ui(x 2 ) is a weak function of X 2 along the vertical span. Of course, ufuj 
is a mirror symmetric reflection of uYuJ. 

Reynolds Stress Closure Verification 

No complete experimental data set exists for quantitative comparison of 
the juncture region turbulent flow 3DPNS prediction. For this geometry, the 
initially steep chordwise pressure gradient is principally responsible for 
generation of the axial vortex pair. Thereafter, in the mid-chord region, 
the combined action of milder pressure gradients and Reynolds stress distribu- 
tions govern the detailed flow-field evolution. The typical experimental juncture 
region geometry, cf. [12, 13], is constructed as the right intersection of two 
plane surfaces with noncoincident leading edge. One surface is the wind tunnel 
floor, while the second is a rounded leading edge, nonfilletted, finite thick- 
ness flat plate mounted perpendicular to the floor. The resultant three-dimen- 
sional, separated stagnation region flow yields a pressure-gradient induced 
roll vortex, in the "wing" leading edge region, which is then convected down- 
stream under nominal zero axial pressure gradient. 

An experimental configuration that specifically facilitates the Reynolds 
stress closure verification is turbulent flow in a straight, uniform rectangular 
cross-section duct, cf. [11]. Following the localized entrance region effects, 
the mild axial pressure gradient is nominally uniform on the cross-section and 
of magnitude sufficient to compensate for duct friction losses. Far downstream, 
experiments verify that no consequential transverse plane velocities exist for 
laminar flow. Conversely, for a turbulent flow, four persistent axial vortex 
pairs exist, one in each right angle corner of the duct. For the experimental 
specification of [11], Baker and Orzechowski [32] document qualitative agree- 
ment for the 3DPNS algorithm prediction on a coarse 13X13 (M = 288) discretiza- 
tion of the symmetric quarter duct. Specifically, no vortex pair roll-up 
occurred for the laminar flew, or for the turbulent flow prediction with the 
0(6] terms involving Ux derivatives in u?uT set to zero, see equation 22. 

Xj jg 
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However, with these terms included, as required theoretically by invariance 
within the constitutive theory, the turbulent flow 3DPNS prediction immediately 
generated the corner vortex pair. Therefore, in this instance of a nonpressure 
driven flow, the anisotropy of the Reynolds stress tensor is principally respon- 
sible for generation of the axial vortex pair. 

With the CMC :3DPNS code now operational on the CDC Cyber/203, refined 
grid and full-duct solutions for the configuration of [11] have been executed 
for quantitative comparison to experimental data [11] for transverse plane 
distributions of u^. and at a downstream station. Figure 9a), from [11], 
is the experimental measurement, of transverse plane velocity distribution 
U£(xi = 37, x^), with u'^/ui= 0.0086. Figure 9b), from [32], is the coarse 
(M = 288) grid 3DPNS solution on the symmetric quarter duct, which exhibits 
essential qualitative agreement with data. However, u*^ /Ui = 0.0010 is a 
factor of eight lower than the data, and large vortex patterns are erroneously 
produced adjacent to both symmetry planes. Figure 9c) is the refined (M = 1052) 
grid solution on the symmetric quarter duct. The qualitatively correct vortex 
patterns nearly fill the section, and u*J/ui = 0.0043 is only a factor of two 
lower than the data. The erroneous vortices remain predicted next to both 
symmetry planes, but their size is substantially reduced in comparison to 
Figure 9b). The combined M = 288 and M = 1!052 solutions confirm that the error 
mechanism causing this local pollution of the solution is a singularity 
in the boundary conditions for the mass conservation harmonic function <|)^(x^), 
where the symmetry plane intersects the wall. Since a velocity component is 
permitted (must occur) parallel to the symmetry plane, but not along the no- 
slip wall, the intersection corresponds to a discontinuous switch from homo- 
geneous Dirichlet to homogeneous Neumann boundary conditions at a corner. 

Removal of this singularity requires the 3DPNS solution domain to span the full 
duct cross-section. The 3DPNS solution executed on a coarse ( M = 1052) grid 
discretization of the entire duct, did predict extinction of the spurious vor- 
tices. However, the extremum transverse velocity = 0.0020 is a factor of 

four lower than data, indicating the discretization too coarse for qualitative 
solution comparison. Hence, a 50 X 50 (M « 2500) discretization is thus 
indicated as the refinement required for a 3DPNS simulation of a ducted 
turbulent flow. 


34 




Figure 9. 3DPMS and Experimental Distributions of Transverse Plane Velocity u^ 

Turbulent Rectangular Duct Flow. 
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) 3DPNS Solution, M = 1052, u'J/Ui= 0.0043 


Figure 9. Concluded 


Reference [11] also documents experimentally measured distributions of 
Ui, k, and u^uj at the downstream station. Figures 10-17 compare these data 
with the 3DPNS solution obtained on the M = 1052 quarter-duct discretization. 
Figure 10a) indicates intrusion of the high momentum core velocity Ui into 
the corner region, as induced by the vortex structure, Figure 9 a). The 3DPNS 
solution exhibits this character for Ui <_ 0.70 only. Figure 10b). The inter- 
section of the 3DPNS Ui = 0.70 isovel with the symmetry plane is in good 
agreement with experiment. Above this level, and on the symmetry planes, the 
intersection of 3DPNS levels for ui exceed data by Ax^^^ = 15%. Along the corner 
bisector, the levels are in better agreement. 

Figures 11-13 compare the Reynolds shear stress distributions. Good over- 
all agreement on level is indicated, as well as some detail of the contour 
shapes for the largest levels. The experimentally measured regions of small 
negative (ITfuJ) and small positive (ujus) shear stress result from the curve 
inflections in ui. Figure 10a), The 3DPNS solution has correctly predicted 
this essential character, although the details are largely affected by the 
boundary condition singularity effects. This is clearly illustrated by the 
"ears" on the 3DPNS prediction of U 2 U 3 , Figure 13. No experimental determin- 
ation of this shear stress component is reported in [ 11 ], since the signal to 
noise ratio at 0 ( 10 ”^*) is essentially unresolvabl e. 

Figures 14-17 summarize the comparison of the square root of the Reynolds 
normal stresses and the turbulence kinetic energy level. Overall, the agree- 
ment on levels of Hi is good, confirming use of the standard definitions for 

the constitutive equation model constants equation 22 , and the coefficients 
6 ^ 

in equations 23-24. The 3DPNS prediction for Uj and k. Figures 14-15, are 

sjmimetric in agreement with data. The intrusion of the lower levels from the 
core region along the bisector is substantially under-predicted. The inter- 
section of TIT = 0.075 on the symmetry plane is in good agreement, and the 
higher 3DPNS solution levels exhibit better agreement with data. The inter- 
section of uT = 0.05 is different by Ax^ x 20%, indicating the level of turbu- 
lence in the experimental core flow is considerably larger than that of the 
3DPNS simulation, see also Figure 15. The 3DPNS solution initialization level 
of k° ~ 10 "® in the potential core region is considerably smaller (probably) 
than the experiment. 
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b) 3DPNS Solution 
Xi/D = 35.8. 


Figure 10. 3DPNS and Experimental Distributions, 
Mean Velocity Ui, Turbulent Rectangular Duct Flow. 
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a) Experimental [11] 
Xi/D = 37.0. 
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b) 3DPNS Solution 
Xi/D = 35.8. 


Figure 11. 3DPNS and Experimental Distributions, Reynolds 
Shear Stress -uTuT X 10® , Turbulent Rectangular Duct Flow. 
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Figure 13. 3DPNS Solution Distribution - Reynolds Shear Stress 
-U 2 U 3 X 10^^, !";ji tiij I pnt Rei .dtiguiar uuct t-iow. 


Figures 16-17 compare the transverse plane normal stress distributions. 
Agreement on overall levels is considerably better, with the 3DPNS prediction 
exhibiting the essential nonsymmetries of the experimental data. Note the 
two 3DPNS solutions are mirror symmetric, while the data are less so. Most 
importantly, recall that these (modest) nonsymmetries are computationally 
confirmed to be the principal causal mechanism of the counter-rotating vortex 
structure for turbulent flow in a straight, rectangular cross-section duct. 
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a) Experimental [11] 
Xi/D = 37.0 . 




Figure 14. 3DPNS and Experimental Distributions, Reynolds 
Normal Stress -/ uYu^, Turbulent Rectangular Duct Flow- 
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Figure 16. 3DPNS and Experimental Distributions, Reynolds 
Normal Stress -/ \ifu{ > Turbulent Rectangular Duct Flow . 
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CONCLUSIONS 


An order of magnitude analysis has yielded a consistent physical variables 
formulation for the parabolic approximation to the three-dimensional Navier- 
Stokes equations for steady, turbulent subsonic flow. A finite element numer- 
ical solution algorithm is derived that accurately enforces the dominant 
differential equation set through formulation of a penalty differential con- 
straint statement. A tensor field expansion is employed to provide closure 
for the Reynolds stress distribution, in concert with solution of two turbulent 
transport equations. A composite pressure field construction is identified to 
enforce overall ellipticity, using a multipass interaction solution procedure 
with a three-dimensional potential flow exterior solution. Numerical results 
document the robustness of the key elements of the developed algorithm for the 
aerodynamic juncture region geometry. 
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APPENDIX 


The finite element algorithm statement for the 3DPNS equation system is 
readily recast into equivalent FORTRA N statements using a hypermatrix formu- 
lational structure [20]. The operation basic to the finite element algorithm 
equation 39 is integration of products of the elements of the cardinal basis 
{N«(x.)}, and the associated (gradient) derivative t:— {Ni(x.)}, for the 

Jo o X - J6 

discretization of R^ formed by the union of triangles. The master element is 
graphed in Figure A.l, which illustrates the various required coordinate 
systems including the linearly dependent natural coordinate system . The 
elements of {Ni> are identical to and for any domain R| 


-P ^r .-y 
^2 ^3 dx 


"e 


Ae p! q! r! 

2 (2 + p + q + r)! 


where A^ is the plane area of R|. Furthermore, 


{?} = 


1 - ni/iii - (1 - ni/rii)f| 2 /ri 2 
< m/fii - (fi!/Ti?)fi2/Til ^ 

m/nl 


(A.l) 


(A. 2) 


and n“ denotes the n„ coordinate of node point a, for the sequencing defined 

A/ 

in Figure A.l. The elementary transformation defining = n^(n|^) is, 




(A. 3) 


where are the direction cosines defining rii‘ as the line connecting nodes 1 
and 2 of R2 . The derivatives of the elements of {Ni(?)> are formed, using 
the chain rule and tensor index summation convention, as 



{Ni(c)} = 


3 


aci 



(A. 4) 
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Consider the first term of the 3DPNS generalized differential equation 
whereby the elements of q-(x. ) = {Uj, U 2 » U 3 , H, k, e} are marched in the 

vl * 

chordwise (x^) direction. For the grid-stretching coordinate transformation, 
equation 31, and subtracting out the continuity equation 1, which yields the 
''nonconservative form", this term becomes 

’^2h23)|^’ + (hz3 + ri3h33)|;^qj (A. 5) 

Employing the finite element construction equation to interpolate pUj and 
nJi(c) on R 2 , yields 

(pui)^ = S{Ni(O}^{RH0Ul}g 
e 

nj = Z{Ni(C)}'''{ETAL} 
e 

The elements of {RH0Ul}g and {ETAL}g are the nodal values of pUi and 
Then, the term corresponding to equation A. 5 within the error extremization 
weighted residuals statement, equation 39, upon rearrangement of selected 
scalars, becomes 




(A. 6 ) 
(A. 7) 


Equation A . 8 defines the global calculus operations on R^ = UR^ , as the 
matrix assembly (S^) of the equivalent calculations performed on the master 
element. (Baker £20, Ch.2] discusses this topic in depth, and rigorously 
derives the matrix row summation procedures on R^ which constitute "assembly") 
Note also that the element matrices {RH0Ul(xi)}g and {QJ(xi)}g are indepen- 
dent of n|^» and can be extracted outside the integral as shown. 
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Considering the first term in equation A. 8, and expanding in terms of {?} 
yields 


{RHOUl}' 


^ {Ni>{Ni}{Ni}^dn {QJ}„ 

"“C 1 1? 2 ^ 3 


= {RHOUl} 


Ki 


^ , C 2^ 3 


dn {QJ}, 


(sym) 


“ 2 _J 


(A. 9) 


Recalling the matrix rules of scalar multiplication, the column matrix 
premultiplier in equation A. 9 can be brought inside by multiplying every 
element of [*J by {•}. This yields a 3 x 3 square matrix, every element of 
which is a 3 x 1 column matrix, i.e., a "hypermatrix" of degree one. The 
defined integrals of products of c^.,1 £ i £ 3, are easy to evaluate using 
equation A.l, which yields the finite element matrix equivalent of equation 
A. 9 as 


(A. 10) 

In equation A. 10, A^ is the element plane area, and 1/60 is the normalizing 
coefficient of the integers constituting the "standard" master hyper-matrix 
[B3000J. This lexicographic symbol indicates the master matrix is defined on 
a two-dimensional element (B), and is constituted of the product of three 
cardinal basis (3), none of which are differentiated (000). If a basis other 
than {Ni(c)} were chosen, the specific entries in [B3000[| would change, of 
course, but the symbolic representation in the final line of equation A. 10 
remains unaltered. Matrix multiplications must clear the hyper-matrix rank 
first; thereafter, the regular rules of matrix algebra apply. 


{RHOUnl f {Ni}{Ni}{Ni}’'’dti{QJ}g= ^{RHOUl}^ 


e 



f6] (r 


[2^ 





1 

► 


W lij 

1 

I 2 J 


f2l 


rii 


^6 

■ 

2 

1. 

I 2 J 


2 

V. > 




(2' 

1 

(sym) 


2 

6 

1 


= A {RH0Ul}g[B3000]{QJ}g 
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Returning. -to equation A. 8, the grid stretc+iing coordinate transformation 
has introduced additional terms which involve derivatives of elements of 
{Ni(?)} on the ti|^. Using equation A. 4, 


|-{Ni(c)} = h2i(xi){B112} §2 + h3i(xi){B113> 63 


(A. 11) 


an 


The elements of {BllK}g, which are element-dependent, 3x1 column matrices, 
are strictly a function of the node coordinates n? of , see equation A. 2, 
and the set a^j^ of direction cosines, equation A. 3. The unit vectors 
are parallel to the x^ coordinate system, and the metrics \l are functions 
of Xj at most. Therefore, equation A. 11 is the matrix equivalent of the 
directional derivative with scalar components parallel to u^ . Then, on the 
master element, the second and third terms in equation A. 8 are of the form 


{RHOUl} 


f{Ni}{Ni}[-hji2- h^3 {Ni}‘^{ETAL}J 




= A^ {RHOUl}' 
e e 


-h22 h2i[B200] - has h2i[B3000]{ETA2}gJ{B112}^{QJ}^ 


+ [-h32 h3i[B200] - h33 h3i[B3000]{ETA2}J{B113}^{QJ}g 


(A. 12) 


The master hyper-matrix [B3000 ]was defined in equation A. 10; it is an ele- 
mentary operation to further show that 


[B200] - 


"2 1 r 
2 1 
(sym) 2 


(A. 13) 
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In distinction to the comments regarding the universality of the symbol 
[B3000] for all (k— degree) cardinal basis specifications the form of 

equation A. 13 would change if the elements of t — { K } are functions of n„- 

dri . K X, 

The general statement for equation A. 8 is, ^ 


{N. Hpul* Ijj) dx = A^{RH0Ul}J[B3000]{QJ}r 

2 K oXi e c 0 

e 


+ {RHOUl} 


[B300L]g - h^3h^i [B400L0]g{ETAL}J{QJ}^ 


(A. 14) 


In equation A. 14, the index £(L) is a tensor index that takes the values 

2 (£,L) £ 3. Furthermore, [B400L0J is a hypermatrix of degree two; the 

first and last Boolean indices (0) indicate {RHOUl} and {ETAL}„ are inter- 

0 6 

polated, and these (inner) multiplications must be performed prior to post- 
multiplication by {QJ}g> which has been differentiated parallel to n^(L). 
The remaining terms in the finite element algorithm statement, equation 39, 
are formed in the same manner 


The second major formulational step is construction of the Jacobian of the 
Newton iteration algorithm, equation 42. By definition. 


r-ii _ 9{FJ} 

" 3{QI} 

with {FJ} given by equation 41. Continuing with the example of the down- 
stream convection term, the specific form for the resultant expressions in 
{FJ} is 

{F.n =s h. , 

e e‘ 


(A. 15) 


A^{RH0U1}^[B3000]({QJ}^^^ - (QJ}j) 


AX i 

2 


Ag{RH0Ul}g 


h£2h£,[B300L] 


h^3hjt,[B400L0] {ETAL} 


{QJ} + 


' j+1 .j 


(A. 16) 


52 


In equation A. 16, {RHOUDg = ^({RH0Ul}P^j +tRH0Ul>.), and AXi = - x^, 

is the chordwise marching step-size. Furthermore, *1 . indicates the 

V+1, J 

algebraic average, and superscript p is the iteration index, see equation 43. 

For equation A. 15, the independent variable is {QJ}P,i • therefore, the elemen- 

J + i 

tal contributions to jj] become formed as 

= A^{RTBU1}^[B3000]«jj 

- ^ {RH0Ul}T£hj,h^i[B300L]^ 

+ lB400L0]{ETAL)1«jj 

, 3{F0> 9{RH0U1} 

alRHOUl} 3{QI} (A. 17) 


In equation A, 17, is the discrete index Kronecker delta, which yields the 
self-dependence expressions, i.e., 9{FJ}/3{QJ}. Since pu^, i.e., {RHOUl} , is 
a function of both p and Ui, which are dependent variables, the second term 
in equation A. 17 yields the nonself coupling. The algebraic equation of 
state yields p = p(p,H); for subsonic flows, the density variation is very 
weak and therefore can be neglected in [J]. Using the chain rule, then 


3pu? 


P'S 11 


(A. 18) 


Hence, interchanging orientation in the hypermatrix formations, as required, 
yiel ds 

airuj diKnuui j _ 

3{QI} " 2 ^"-'e‘ 


{QI}I[B3000]5jj6j3^ 


3 {RHOUl} 


AXiAgPe 


{QI>l|h^2h^i[B3000] 


+ 


h£3h£i[B4L000]{ETAL}g 


'^Jl'^Il 


(A. 19) 
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In equation A. 19, is the element-average value of p on R|. The CMC:3DPNS 
computer code neglects all nonself coupling in construction of the Jacobian 

i- 

|J]. Therefore, [J][ is independent of the specific element of q^. = {Ui, U 2 . 

U3, H, k,e}. Hence, the single LU decomposition of [J(QJ)]/+i ''S employed to 

D"l" 1 0 

solve for the appropriate six elements of using a multiple right-hand 

side procedure in equation 42, The Jacobian is updated for each iteration 
within step 

Equations A. 14, A. 17 and A. 19 are illustrative of the operational procedures 
of formulation of the finite element algorithm statement for aerodynamic 
juncture region flow. There is an exacting amount of detail required to 
complete all aspects of the 3DPNS algorithm statement. However, the developed 
hypermatrix formalisms and master element concepts have produced a rigorous 
procedure to keep track of the details. The role of the tensor indices and 
matrix differential calculus are invaluable tools put to practical use. 

Finally, note that these equations are written in a pseudo-FORTRAN language 
which yields coding which is nominally identical in appearance to the theo- 
retical statements. 
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